#define CUDA
#ifdef CUDA
#include <hip/hip_runtime.h>
#include <hipblas.h>
#endif


#include "norm.hpp"
#include "gcr.h"
#include <math.h>
#include <iostream>
#include <fstream>
#include <mpi.h>
//#include <mkl.h>
#include <mkl_cblas.h>
#include "gptl.h"
#include "matrixpro.hpp"
#include "svrasr.hpp"
#define new_max 20
#define precondition
#define mpi_reduce
#define index_h(i,j) ((i) * (iter_max) + (j))
#define index_g(i,j) ((i) * NG1 + (j))
#define index_v(i,j) ((i) + (j) * (iter_max+1))
#define index_kr(i,j) ((i)+(j)*(kR+1))
//#define index_h(i,j) (i +(j)*iter_max)
extern "C" void matrixpro_c(float *a,float *b,float *c,int its, int ite, int jts, int jte, int kts, int kte, int jend);
extern "C" void matrixpro_x(float *a,double *b,float *c,int its, int ite, int jts, int jte, int kts, int kte, int jend, int,int,int,int);
extern "C" void module_halo_mp_glob_updatehalo_real_3d_(float* h_p, int *,int *); 
extern "C" void svrasr( float*  yl,int iter_max, float*  b, int ni, int nk, int nj);
extern "C" void glob_updatehalo_p_(float *h_p, int *m, int *iter_max,int *its, int *ite, int *kts, int *kte, int *jts, int *jte);
#define printH()\
    fprintf(fp,"\nH=\n");\
        for(int gi=0;gi<=iter_max;gi++){\
            for(int gj=gi;gj<iter_max;gj++){\
                fprintf(fp,"%3.2E  ",H[index_h(gi,gj)]);\
            }\
            fprintf(fp,"\n");\
        }
//#define printH0()\
    fprintf(fp,"\nH0=\n");\
        for(int gi=0;gi<=new_max;gi++){\
            for(int gj=0;gj<new_max;gj++){\
                fprintf(fp,"%5.4E  ",H0[index_h(gi,gj)]);\
            }\
            fprintf(fp,"\n");\
        }
//#define printH()\
    if(tid==0){\
    printf("\nH=\n");\
        for(int gi=0;gi<=new_max;gi++){\
            for(int gj=0;gj<new_max;gj++){\
                printf("%5.4E  ",H[index_h(gi,gj)]);\
            }\
            printf("\n");\
        }}

#define printH0()\
    if(tid==0){\
    printf("\nH0=\n");\
        for(int gi=1;gi<=new_max;gi++){\
            for(int gj=gi-1;gj<new_max;gj++){\
                printf("%3.2E  ",H0[index_h(gi,gj)]);\
            }\
            printf("\n");\
        }}
#define printQ()\
if(tid==0){\
    printf("\nline%d,Q=",__LINE__);\
        for(int gi=0;gi<=new_max;gi++){\
            for(int gj=0;gj<new_max;gj++){\
                printf("%4.2E  ",Q[index_kr(gi,gj)]);\
            }\
            printf("\n");\
        }\
}
#define printQ1()\
if(tid==0){\
    printf("\nQ=");\
        for(int gi=0;gi<=kR;gi++){\
            for(int gj=0;gj<=kR;gj++){\
                printf("%4.3E  ",Q1[index_kr(gi,gj)]);\
            }\
            printf("\n");\
        }\
}


#define printLOO()\
                c1[0]=cblas_sdot(iter_max*iter_max,Q1,1,Q1,1);\
                MPI_Allreduce(c1,c2,1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);\
                if(tid==12) printf("mytid=%d,LOO=%E,tol=%E\n",tid,sqrt(c1[0]),sqrt(c2[0]));
double norm2(const float *v, const  int n)
{
    double tmp = 0;
    #pragma omp parallel for  reduction(+:tmp) schedule(static)
    for(int i=0; i<n; i++)                                                                                                                
        tmp += v[i] * v[i];
    return tmp;//sqrt(tmp);
}
__global__ void init_p2(float *d_p, float *d_r,  int its, int ite, int jts, int jte, int kts, int kte,int jend) {
    #include "defindex.hpp"
    cupara_jte(d_p[index4b(i,k,j,0)]=d_r[index3(i,k,j)];)
}
__global__ void init_p1(float *d_p,  float *d_r,  float beta, int its, int ite, int jts, int jte, int kts, int kte,int jend) {
    #include "defindex.hpp"
    cupara_jte(d_p[index4(i,k,j,0)]=d_r[index3(i,k,j)]*beta;)
    //cupara_jte(d_r0[index3b(i,k,j)]=d_r[index3(i,k,j)];)
}
__global__ void init_p(float *d_p, double *d_x0, int its, int ite, int jts, int jte, int kts, int kte,int jend) {
    #include "defindex.hpp"
    cupara_jte1(d_p[index4b(i,k,j,0)]=d_x0[index_x(i,k,j)];)
}
__global__ void trans_a(float* __restrict__ d_a_helm, float*  __restrict__ d_a0, int its, int ite, int jts, int jte, int kts, int kte,int jend) {
    #include "defindex.hpp"
    for(int m=1;m<=19;m++){
        cupara_jte(d_a0[index_a(m,i,k,j)]=d_a_helm[index_a0(m,i,k,j)];)
    }
}
__global__ void gcr_init( float* __restrict__ d_f0, float* __restrict__ d_r,int its, int ite, int jts,int jte,int kts, int kte, int jend){
	#include "defindex.hpp"
        cupara_jend(d_r[index3(i,k,j)]  = d_f0[index3(i,k,j)]-d_r[index3(i,k,j)];)
//        cupara_jend(d_r0[index4b(i,k,j,0)] =d_r[index3(i,k,j)];)
	}

__global__ void gcr_iteration(int iter_max,float* __restrict__ d_a0, float* __restrict__ d_p, float* __restrict__ d_ap,int its, int ite, int jts,int jte,int kts, int kte, int jend){
	#include "defindex.hpp"
        cupara_jte1(d_p[index4b(i,k,j,0)] = d_p[index4b(i,k,j,iter_max-1)];)
}
__global__ void gcr_iteration4(int step,  float* __restrict__ d_xc,  float* __restrict__ d_p,  float* __restrict__ d_ap, float* __restrict__ d_p1,  int its, int ite, int jts,int jte,int kts, int kte, int jend){  
    #include "defindex.hpp"
    cupara_jend(
//        d_p1[index4b(i,k,j,0)]= d_xc[step]* d_p[index4b(i,k,j,step)];
        for(int l=0;l<=step;l++){
            //d_p1[index4b(i,k,j,0)]-=d_xc[l]*d_p[index4(i,k,j,l)];
            d_ap[index3(i,k,j)]-=d_xc[l]*d_p[index4(i,k,j,l)];
        }
        //d_x0[index_x(i,k,j)] += d_p1[index4b(i,k,j,0)];
        )
}


__global__ void gcr_iteration4b(int step, float beta, float* __restrict__ d_xc,  float* __restrict__ d_p, float* __restrict__ d_p1,  int its, int ite, int jts,int jte,int kts, int kte, int jend){  
    #include "defindex.hpp"
    cupara_jte1(
            d_p1[index4b(i,k,j,0)]=0.0;
            )
    cupara_jte(
//        d_p1[index4b(i,k,j,0)]= d_xc[step]* d_p[index4b(i,k,j,step)];
        for(int l=0;l<=step;l++){
            d_p1[index4b(i,k,j,0)]+=d_xc[l]*d_p[index4(i,k,j,l)];
        }
        //d_x0[index_x(i,k,j)] += d_p1[index4b(i,k,j,0)];
        )
//    cupara_jte(d_p1[index4b(i,k,j,0)] *= beta;)
}
__global__ void gcr_iteration3b(int iter_max, int kR, float* __restrict__ d_p, float* __restrict__ d_xc,float* __restrict__ d_P, int its, int ite, int jts,int jte,int kts, int kte, int jend){
		#include "defindex.hpp"
            //initialize(G,NG1*(kR+1));
            //for(int l=0;l<=iter_max;l++){
            for(int t=0;t<kR;t++){
                cupara_jte(d_p[index4(i,k,j,kR)]-=d_xc[t]*d_p[index4(i,k,j,t)];)
            }
            //}
}


//hipLaunchKernelGGL(gcr_iteration, numBlocks, threadsPerBlock, 0, 0, step, d_x0, d_xc, d_ap,  d_p, d_p1 , its, ite, jts, jte, kts, kte, jend);
__global__ void gcr_iteration3(int iter_max, int kR, float* __restrict__ d_G,float* __restrict__ d_p, float* __restrict__ d_P, int its, int ite, int jts,int jte,int kts, int kte, int jend){
		#include "defindex.hpp"
            //initialize(G,NG1*(kR+1));
            for(int l=0;l<=iter_max;l++){
                for(int t=0;t<=kR;t++){
                    cupara_jte(d_G[index4(i,k,j,t)] += d_p[index4(i,k,j,l)] * d_P[index_v(l,t)] ;)
                }
            }
/*
        for(int l=0;l<=m;l++){
            double bl=-c2[l]/aps[l];
            cupara_jte1( 
                    d_p[index4b(i,k,j,m+1)]=d_p[index4b(i,k,j,m+1)]+bl*d_p[index4b(i,k,j,l)];
                    )
            } */
}
/*
__global__ void gcr_iteration2(int m,float* __restrict__ d_a0,double *d_x0, float* __restrict__ d_ar,double  d_ac, float* __restrict__ d_p, float* __restrict__ d_r, float* __restrict__ d_ap,int its, int ite, int jts,int jte,int kts, int kte, int jend){
	#include "defindex.hpp"
    cupara_jend(
        d_x0[index_x(i,k,j)]= d_x0[index_x(i,k,j)]+d_ac*d_p[index4(i,k,j,m)];
        d_r[index3(i,k,j)] = d_r[index3(i,k,j)]-d_ac*d_ap[index4(i,k,j,m)];
    )
    cupara_jte(d_p[index4b(i,k,j,m+1)] =d_r[index3(i,k,j)];)
}
*/
__global__ void gcr_final(double* __restrict__ d_x0, float* __restrict__ d_ap, float* __restrict__ d_r,  float* __restrict__ d_p1, int its, int ite, int jts,int jte,int kts, int kte, int jend){
    #include "defindex.hpp"
        //float beta1=1/beta;
        //cupara_jend(d_r[index3(i,k,j)]  *= beta;)
        cupara_jend(d_r[index3(i,k,j)]  -= d_ap[index4(i,k,j,0)];)
        cupara_jte1(d_x0[index_x(i,k,j)] += d_p1[index4b(i,k,j,0)];)
        }

/*__global__ void gcr_final1(double* __restrict__ d_x0, float* __restrict__ d_r, float* __restrict__ d_f0, float* __restrict__ d_p1, int its, int ite, int jts,int jte,int kts, int kte, int jend){
    #include "defindex.hpp"
        cupara_jend(d_r[index3(i,k,j)]  = d_f0[index3(i,k,j)] - d_r[index4(i,k,j,0)];)
        cupara_jte1(d_x0[index_x(i,k,j)] += d_p1[index4b(i,k,j,0)];)
        }*/
double max_norm(const float *v, const  int n)
{
    double tmp = 0, max_val=0;
    for(int i=0; i<n; i++){
        tmp=fabs(v[i]);
        if(tmp>max_val) max_val=tmp;
    }
    return max_val;//sqrt(tmp);
}




//! original update operations
void Update(float *y, const int k,  float *H, const int iter_max, float *rhs)
{
    int i=0, j=0;
    //float *y = (float*) malloc((k+1)*sizeof(float));
    for (i=0; i<=k; i++)  
        y[i] = rhs[i];

    // Back substituation, H is upper triangle matrix
    // solve y, where H*y = rhs
    for (i = k; i >= 0; i--) {
        for (j=i+1; j <= k; j++)
            y[i] -= H[index_h(i,j)] * y[j];
        y[i] /= H[index_h(i,i)];// the diagonal element
    }
}

extern "C" void mygivens(float a0, float b0, float* sn, float* cs, double &r, int i) {
    float a=a0, b=b0;
    r = sqrt(a*a + b*b);
    cs[i] = a/r;
    sn[i] = b/r;
}


//Q = B*inv(A) where A is upper triangular.
void rt_mult_inv(float *B, float *A, float *Q, int js, int iter_max){
    
    for(unsigned int i = 0;i<=js+1;i++){//col
        Q[index_h(i,0)] = B[index_h(i,0)]/A[index_h(0,0)];
        for(unsigned int j = 1; j<=js;j++){
            double temp = 0; 
            for(unsigned int l = 0; l<j;l++){
                temp += Q[index_h(i,l)]*A[index_h(l,j)];
            }    
            Q[index_h(i,j)] = (B[index_h(i,j)]-temp)/A[index_h(j,j)];
        }    
    }                                                                           
}

void com_LOO(float *Q, float *Q1, int dim, int kR){
        int i,j,k;
    // Compute I-Q^T*Q
    for (i = 0; i <= dim; i++) {
        for (j = 0; j <= dim; j++) {
            double sum = 0;
            for (k = 0; k <= dim; k++) {
                sum += Q[index_kr(k, i)] * Q[index_kr(k, j)];
            }
            Q1[index_kr(i, j)] = i == j ? 1 : 0;
            Q1[index_kr(i, j)] -= sum;
        }
    }

}

extern "C" void cu_glob_updatehalo_p(float* h_p, float* d_p, float* h_buff, float* d_buff,
    int its, int ite, int kts, int kte, int jts, int jte, 
        int m, int iter_max, int* mytid);


//#define share_mem

//void gmres_iter(float *a0, float *H, float *H0, float *h_p, float *h_zp, float *h_ap, float *cm_helm, float *rhs,float *rhsR, int &js, int iter_max, int its, int ite, int jts, int jte, int kts, int kte, int jend, int tid, int kR, float *Q, float ep, int &flag, float *h_r, float *h_r0){
void gmres_iter(float *a0,float* d_a0, float *H, float *H0, float*d_p, float *h_zp, float*d_p1,  float*d_ap, float*d_ap1,   float *cm_helm, float*d_cm_helm, float *rhs,float *rhsR, int &js, int iter_max, int its, int ite, int jts, int jte, int kts, int kte, int jend, int tid, int kR, float *Q, float ep, int &flag, float*d_r, float * buff,float * d_buff, float *d_xc, hipblasHandle_t h){
        int NX=ite-its+3, NY=kte-kts+3, NZ=jte-jts+3, NG1=NX*NY*NZ, NG=(NX-2)*NY*(NZ-2);
        //float c1[iter_max+1]; 
    dim3 threadsPerBlock(64, 4, 1);
    dim3 numBlocks(ceil(((float)NX)/64),ceil(((float) NY)/4),NZ);
    dim3 threadsPerBlock2(64, 1, 1); 
    dim3 numBlocks2(1,NZ,1); //halo=1

    dim3 block_size(64,4,1);
    dim3 num_blocks((NX + block_size.x - 1) / block_size.x,(NY + block_size.y - 1) / block_size.y,NZ);
        #ifdef share_mem
        dim3 threadsPerBlock_mat(16, 4, 6);
        dim3 numBlocks_mat(ceil(((float)NX)/(threadsPerBlock_mat.x-2)),ceil(((float) NY-1)/(threadsPerBlock_mat.y)),ceil(((float)NZ)/(threadsPerBlock_mat.z-2)));
        #else
        dim3 threadsPerBlock_mat(16, 4, 4);
        dim3 numBlocks_mat(ceil(((float)NX)/16),ceil(((float) NY)/4),ceil(((float)NZ)/4));
        #endif
//    FILE* fp = fopen("gmres.txt", "w");
    float *c1 =new float[iter_max+1];
    float *c2 =new float[iter_max+1];
    int i,j,k, is, zero=0;
    double tmp_sum=0.0;
    float *cs =new float[iter_max+1];
    float *sn =new float[iter_max+1];
    initialize(cs,iter_max+1);
    initialize(sn,iter_max+1);
    double relres;
    for(js=kR;js<iter_max;js++){
        GPTLstart ("gmres_iter");  
        //CHECK(hipMemcpy(d_p+js*NG1,h_p+js*NG1,NG1*sizeof(float),hipMemcpyHostToDevice));
    hipLaunchKernelGGL(init_p2, num_blocks, block_size, 0, 0, d_p1,d_p+js*NG, its, ite, jts, jte, kts, kte, jend);

        cu_glob_updatehalo_p(h_zp, d_p1, buff, d_buff,its, ite, kts, kte, jts, jte, js, iter_max,  &tid) ;
        GPTLstart ("precondition");  
        hipLaunchKernelGGL(pacu_svrasr, numBlocks2, threadsPerBlock2, 0, 0, d_p1, 0,d_cm_helm, NX, NY, NZ);
        cu_glob_updatehalo_p(h_zp, d_p1, buff, d_buff,its, ite, kts, kte, jts, jte, zero, iter_max,  &tid) ;
        GPTLstop ("precondition");  

        GPTLstart ("matrixpro");  
        //CHECK(hipDeviceSynchronize());
        hipLaunchKernelGGL(cu_matrixpro, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_p1,d_ap, its, ite, jts, jte, kts, kte, jend);
        CHECK(hipDeviceSynchronize());
        GPTLstop ("matrixpro");  
        GPTLstart ("dot_p");  
        //CHECK(hipDeviceSynchronize());
        for (is = 0; is <= js; is++) {
            hipblasSdot(h, NG,d_p+is*NG,1, d_ap,1, d_xc+is);
        }
            CHECK(hipMemcpy(c1,d_xc, (js+1)*sizeof(float), hipMemcpyDeviceToHost));
        //CHECK(hipDeviceSynchronize());
        GPTLstop ("dot_p");  
            #ifdef mpi_reduce
            MPI_Allreduce(c1,c2,js+1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
            #endif
        GPTLstart ("update_p");  

        CHECK(hipMemcpy(d_xc,c2, (js+1)*sizeof(float),hipMemcpyHostToDevice));
        hipLaunchKernelGGL(gcr_iteration4, numBlocks, threadsPerBlock, 0, 0, js,  d_xc,  d_p, d_ap, d_ap1 , its, ite, jts, jte, kts, kte, jend);
        for (is = 0; is <= js; is++) {
            H[index_h(is,js)] = c2[is];  
        } 
    	//CHECK(hipDeviceSynchronize()); //test time
        GPTLstop ("update_p");  
#define origin
#ifdef origin
        GPTLstart ("normalize_p");  
        hipblasSdot(h,NG,d_ap,1,d_ap,1,d_xc);
        CHECK(hipMemcpy(c1,d_xc, 1*sizeof(float), hipMemcpyDeviceToHost)); 
        MPI_Allreduce(c1,c2,1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
        tmp_sum=sqrt(c2[0]); H[index_h(js+1,js)]=tmp_sum;
        //CHECK(hipMemcpy(d_p+(js+1)*NG,h_ap,NG*sizeof(float),hipMemcpyHostToDevice));
        hipLaunchKernelGGL(init_p1, num_blocks, block_size, 0, 0, d_p+(js+1)*NG, d_ap,(float)1.0/tmp_sum, its, ite, jts, jte, kts, kte, jend);
/*
        hipblasSdot(h,NG,d_r,1,d_ap,1, d_xc+1);
        CHECK(hipMemcpy(c1,d_xc, 2*sizeof(float), hipMemcpyDeviceToHost)); 
        MPI_Allreduce(c1,c2,2,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
*/
        c2[1]=0.0;
        rhs[js+1]=c2[1]/tmp_sum; rhsR[js+1]=rhs[js+1];

#else

        GPTLstart ("normalize_p");  
        hipblasSdot(h,NG,d_ap,1,d_ap,1,d_xc);
        hipblasSdot(h,NG,d_r,1,d_ap,1, d_xc+1);
        CHECK(hipMemcpy(c1,d_xc, 2*sizeof(float), hipMemcpyDeviceToHost)); 
        MPI_Allreduce(c1,c2,2,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
        tmp_sum=sqrt(c2[0]); H[index_h(js+1,js)]=tmp_sum;
        hipLaunchKernelGGL(init_p1, num_blocks, block_size, 0, 0, d_p+(js+1)*NG, d_ap,(float)1.0/tmp_sum, its, ite, jts, jte, kts, kte, jend);
        //CHECK(hipMemcpy(d_p+(js+1)*NG,h_ap,NG*sizeof(float),hipMemcpyHostToDevice));
        //c2[1]=0.0;
        rhs[js+1]=c2[1]/tmp_sum; rhsR[js+1]=rhs[js+1];
#endif
        GPTLstop ("normalize_p");  
        for(is=0;is<=js+1;is++)H0[index_h(is,js)]=H[index_h(is,js)];
            
        if(kR!=0){
            for(is=0;is<=kR;is++){
                tmp_sum=0.0;
                for(int l=0;l<=kR;l++){
                    tmp_sum+=Q[index_kr(l,is)]*H0[index_h(l,js)];
                }
                H[index_h(is,js)]=tmp_sum;
            }
        }
        for (is = kR; is < js; is++){
            tmp_sum=cs[is]*H[index_h(is,js)]+sn[is]*H[index_h(is+1,js)];
            H[index_h(is+1,js)]=-sn[is]*H[index_h(is,js)]+cs[is]*H[index_h(is+1,js)];
            H[index_h(is,js)]=tmp_sum;
        }
        mygivens( H[index_h(js,js)], H[index_h(js+1,js)], sn, cs, tmp_sum,js); 
        H[index_h(js,js)] = tmp_sum;
        H[index_h(js+1,js)] = 0;

        tmp_sum=cs[js]*rhs[js]; rhs[js+1]=-sn[js]*rhs[js];
        rhs[js]=tmp_sum;
        
        relres = fabs(rhs[js+1]);/// normb;
        
        GPTLstop ("gmres_iter");  
        //if(tid==12)
        //    printf("iter %d RES of gmres is %e \n",js,relres);
        
        if (relres  < ep ) {
            if(tid==12)
            std::cout<<"HOST--"<<tid<<"-relres: "<<relres<<",fabs(rhs[i+1])="<<fabs(rhs[js+1])<<std::endl;
            js++;
            flag=0;
            break;
        }
    }
    js--;
    delete[] c1;
    delete[] c2;
    //fclose(fp);
}


extern "C" void qr(float *A, int M, int N, int iter_max, float *Q, float * R);

extern "C" void orth_mgs(float *Vr, int m, int n, char *dreorth, float *P, int tid); 
extern "C" void orth_vec(float *P, float *rR, int m, int n,  int tid); 
//orth_vec(P, rR, iter_max,  kR, tid); 

void swapArrays(float* &a, float* &b) {
    float* temp = a;  // 创建临时指针，并将其指向a指向的数组地址
    a = b;          // 将a指向b指向的数组地址
    b = temp;       // 将b指向临时指针所指向的数组地址
}

extern "C" void  harm_Ritz_vec( int kR, float *H, int n, int *ord, float *vr, int flag, int tid, float *rR); 

#define gmres
extern "C" void  cu_psolve_gmres_main(float ep, float *a0, float *f0, float *cm_helm,int iter_max, double *x0, 
        int dep, int  jdep, int  ids, int  ide, int  jds, int  jde, int  kds, int  kde, 
        int ims, int  ime, int  jms, int  jme, int  kms, int  kme, 
        int its, int  ite, int  jts, int  jte, int  kts, int  kte, int *mytid, int initflag){
    // print_line()
    iter_max=20;
    //const int kR0=4; //4->0
    int L=0;
    int kR=0;
    ep=sqrt(ep); //ep=1e-6; 
    if(mytid[0]==0) printf("\nkR=%d, L=%d, ep=%.3e\n",kR, L, ep);
    int tid=mytid[0];
    //if(mytid[0]==0)printf("ep=%E\n",ep);
#include "defvar.hpp"
        FILE *fp;                                  
        char filename[11];

        sprintf(filename, "file_%d.txt", iter_max);
    //    fp = fopen(filename, "w");     
    //glob_updatehalo_x0_(x0, &ims, &ime, &jms, &jme, &kts, &kte);

    int zero=0,flag=1,is=0,js=0,n=NG;

    //double *y=new double[NG1];
    int *ord=new int[iter_max+1];
    static float *rR =new float[iter_max+1];
    static float *Vr =new float[(iter_max+1) * kR];
    static float *P =new float[(iter_max+1) * (kR+1)];
    float *G =new float[NG1 * (kR+1)];

    float *vr =new float[iter_max * iter_max];
    float *rhs =new float[iter_max+1];// (float*) malloc((iter_max+1)*sizeof(float));
    float *rhsR =new float[iter_max+1];// (float*) malloc((iter_max+1)*sizeof(float));
    initialize(vr,iter_max * iter_max);
    initialize(rhs,iter_max+1);
    static float *H =new float[(iter_max+1)*iter_max];// (float*) malloc(((iter_max+1)*iter_max)*sizeof(float))
    static float *H0 =new float[(iter_max+1)*iter_max];// (float*) malloc(((iter_max+1)*iter_max)*sizeof(float));
    
    if(initflag==1){
    printH0();
//    return;
    }else{
        initialize(H,(iter_max+1)*iter_max);
        initialize(H0,(iter_max+1)*iter_max);
    }
    float *x_c =new float[iter_max];// (float*) malloc((iter_max+1)*sizeof(float));
    static float *Q=new float[(kR+1)*(kR+1)];
    float *Q1 =new float[iter_max*(iter_max+1)];
    initialize(Q1,(iter_max+1)*iter_max);
    float *R=new float[kR*(kR+1)];
    //float normb = norm2(f0, n);
    double relres,beta;
    int outer,maxit=540/iter_max;
    maxit=20; 
    double tmp_sum=0.0; float* c21=new float[kR+1];
    dim3 threadsPerBlock(64, 4, 1);
    dim3 numBlocks(ceil(((float)NX)/64),ceil(((float) NY)/4),NZ);
    dim3 threadsPerBlock2(64, 1, 1); 
    dim3 numBlocks2(1,NZ,1); //halo=1

    dim3 block_size(64,4,1);
    dim3 num_blocks((NX + block_size.x - 1) / block_size.x,(NY + block_size.y - 1) / block_size.y,NZ);
	CHECK(hipMemcpy(d_x0,x0,NG2*sizeof(double),hipMemcpyHostToDevice));
    if(initflag==0){
	CHECK(hipMemset(d_p,0,NG1*(iter_max+1)*sizeof(float)));
	CHECK(hipMemset(d_p1,0,NG1*sizeof(float)));
	CHECK(hipMemset(d_ap1,0,NG1*sizeof(float)));
#ifdef trans
    GPTLstart ("trans_a");
    	//CHECK(hipDeviceSynchronize()); //test time
        hipLaunchKernelGGL(trans_a, num_blocks, block_size, 0, 0, d_a_helm, d_a0, its, ite, jts, jte, kts, kte, jde);
    	CHECK(hipDeviceSynchronize()); //test time
        //initialize(h_p,NG1);
    GPTLstop ("trans_a");
#endif

    //hipLaunchKernelGGL(init_p, num_blocks, block_size, 0, 0, d_p1, d_x0, its, ite, jts, jte, kts, kte, jend);
    hipLaunchKernelGGL(cu_matrixpro_x, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_x0,d_r, its, ite, jts, jte, kts, kte, jend, ims, ime, jms, jme);
    hipLaunchKernelGGL(gcr_init, numBlocks, threadsPerBlock, 0, 0, d_f0, d_r, its, ite, jts, jte, kts, kte, jend);
    hipblasSdot(h, NG,d_r,1,d_r,1, d_xc);                                                                                                            
    CHECK(hipMemcpy(c1, d_xc, sizeof(float), hipMemcpyDeviceToHost));
    beta= c1[0];

    //beta = norm2(h_r, n);
    #ifdef mpi_reduce
    MPI_Allreduce(&beta,&tmp_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    beta=sqrt(tmp_sum);
    #else
    beta=sqrt(beta);
    #endif
    //initialize(h_p,NG1*(iter_max+1)); initialize(h_r0,NG1);
    cuinitialize(d_p,NG1*(iter_max+1)); //initialize(h_r0,NG1);
    hipLaunchKernelGGL(init_p1, num_blocks, block_size, 0, 0, d_p, d_r, (float)1.0/beta, its, ite, jts, jte, kts, kte, jend);
    //hipLaunchKernelGGL(init_p2, num_blocks, block_size, 0, 0, d_r0, d_p, its, ite, jts, jte, kts, kte, jend);
    initialize(rhs,iter_max+1);initialize(rhsR,iter_max+1); rhs[0] = beta;rhsR[0] = beta;
    
#ifdef gmres
    gmres_iter(a0,d_a0, H, H0, d_p, h_zp, d_p1, d_ap, d_ap1, cm_helm, d_cm_helm, rhs, rhsR, js, iter_max, its, ite, jts, jte, kts, kte, jend, tid, 0, NULL, ep, flag,  d_r,  buff, d_buff, d_xc, h );
#endif


    relres = fabs(rhs[js+1]);/// normb;


    Update(x_c, js, H, iter_max, rhs);
    for(int l=0;l<=iter_max;l++){ 
        tmp_sum=0.0;
        for(i=0;i<iter_max;i++)
            tmp_sum+=H0[index_h(l,i)]*x_c[i];
        rR[l]=rhsR[l]-tmp_sum;
    } 
    //initialize(h_zp,NG1);
    cuinitialize(d_p1,NG1);
    CHECK(hipMemcpy(d_xc,x_c, iter_max*sizeof(float),hipMemcpyHostToDevice));
    hipLaunchKernelGGL(gcr_iteration4b, numBlocks, threadsPerBlock, 0, 0, js, 1, d_xc,  d_p, d_p1 , its, ite, jts, jte, kts, kte, jend);

    cu_glob_updatehalo_p(h_zp, d_p1, buff, d_buff,its, ite, kts, kte, jts, jte, zero, iter_max,  mytid) ;
    hipLaunchKernelGGL(pacu_svrasr, numBlocks2, threadsPerBlock2, 0, 0, d_p1, 0,d_cm_helm, NX, NY, NZ);
    cu_glob_updatehalo_p(h_zp, d_p1, buff, d_buff,its, ite, kts, kte, jts, jte, zero, iter_max,  mytid) ;
    hipLaunchKernelGGL(cu_matrixpro, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_p1,d_ap, its, ite, jts, jte, kts, kte, jend);
    hipLaunchKernelGGL(gcr_final, numBlocks, threadsPerBlock, 0, 0,  d_x0, d_ap, d_r, d_p1, its, ite, jts, jte, kts, kte, jend); 
    harm_Ritz_vec(kR,H0, iter_max, ord, vr, 0, tid, rR); // get harmonic ritz vector $vr$
    //beta = norm2(h_r, n);
    hipblasSdot(h, NG,d_r,1,d_r,1, d_xc);
    hipblasSdot(h, NG,d_p1,1,d_p1,1, d_xc+1);
    CHECK(hipMemcpy(c1, d_xc, 2*sizeof(float), hipMemcpyDeviceToHost));
    //beta=c1[0]; //error
    MPI_Allreduce(c1,c2,2,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
    c2[0]=sqrt(c2[0]);
    c2[1]=sqrt(c2[1]);

        if(tid==5) printf("RES of gmres is %e,%e \n",c2[0],c2[1]);
        /*
        fprintf(fp,"\nrhs=\n");
        for(int gi=0;gi<=iter_max;gi++){
            fprintf(fp,"%e ",rhs[gi]);
        }
        fprintf(fp,"\n");
        fprintf(fp,"\nrhsR=\n");
        for(int gi=0;gi<=iter_max;gi++){
            fprintf(fp,"%e ",rhsR[gi]);
        }
        fprintf(fp,"\n");*/
//return;
    for(j=0;j<kR;j++){
        for(i=0;i<iter_max;i++){
            Vr[index_v(i,j)]=vr[index_h(i,ord[j])];
        }
        Vr[index_v(iter_max,j)]=0;
    
    } 
    }else{
        /*
        hipLaunchKernelGGL(init_p, num_blocks, block_size, 0, 0, d_p1, d_x0, its, ite, jts, jte, kts, kte, jend);
        hipLaunchKernelGGL(cu_matrixpro, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_p1,d_r, its, ite, jts, jte, kts, kte, jend);
        hipLaunchKernelGGL(gcr_init, numBlocks, threadsPerBlock, 0, 0, d_f0, d_r, its, ite, jts, jte, kts, kte, jend);
        
        for (is = 0; is <= iter_max; is++) {
        hipblasSdot(h,NG,d_r,1,d_p+is*NG,1, d_xc+1);
        }
        CHECK(hipMemcpy(c1,d_xc, (iter_max+1)*sizeof(float), hipMemcpyDeviceToHost)); 
        MPI_Allreduce(c1,rhsR,iter_max+1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);

        //printH();
        for(is=0;is<(iter_max+1)*iter_max;is++)H[is]=H0[is];
        if(kR>0){
            for(is=0;is<=kR;is++){
                tmp_sum=0.0;
                for(js=0;js<=kR;js++){
                    tmp_sum+=Q[index_kr(js,is)]*rhsR[js];
                }
                rhs[is]=tmp_sum;
            }

            for(js=0;js<iter_max;js++){
            for(is=0;is<=kR;is++){
                tmp_sum=0.0;
                for(int l=0;l<=kR;l++){
                    tmp_sum+=Q[index_kr(l,is)]*H0[index_h(l,js)];
                }
                H[index_h(is,js)]=tmp_sum;
            }
            }
        }
        //printH();
        //return;
        float sn,cs;
        for (is = kR; is < iter_max; is++) {
            //mygivens( H[index_h(is,is)], H[index_h(is+1,is)], sn, cs, tmp_sum,is); 
            tmp_sum = sqrt(H[index_h(is,is)]*H[index_h(is,is)] + H[index_h(is+1,is)]*H[index_h(is+1,is)]);
            cs = H[index_h(is,is)]/tmp_sum;
            sn = H[index_h(is+1,is)]/tmp_sum;

            for (js = is; js < iter_max; js++){
                tmp_sum=cs*H[index_h(is,js)]+sn*H[index_h(is+1,js)];
                H[index_h(is+1,js)]=-sn*H[index_h(is,js)]+cs*H[index_h(is+1,js)];
                H[index_h(is,js)]=tmp_sum;
                tmp_sum=cs*rhs[js]; rhs[js+1]=-sn*rhs[js]; rhs[js]=tmp_sum;
            }
        }
            for(js=0;js<iter_max;js++){
            for(is=js+1;is<=iter_max;is++){
                H[index_h(is,js)]=0.0;
            }}
        //printH();
        Update(x_c, iter_max-1, H, iter_max, rhs);
        
        if(kR>0){
            harm_Ritz_vec(kR,H0, iter_max, ord, vr,1, tid, rR); 
          //
            fprintf(fp,"\n628rR=\n");
            for(int gi=0;gi<=iter_max;gi++){
                fprintf(fp,"%e ",rR[gi]);
            }
            for(int l=0;l<=iter_max;l++){ 
                tmp_sum=0.0;
                for(i=0;i<iter_max;i++)
                    tmp_sum+=H0[index_h(l,i)]*x_c[i];
                rR[l]=rhsR[l]-tmp_sum;
            }
            fprintf(fp,"\n633rR=\n");
            for(int gi=0;gi<=iter_max;gi++){
                fprintf(fp,"%e ",rR[gi]);
            }

            
            for(j=0;j<kR;j++){
                for(i=0;i<iter_max;i++){
                    Vr[index_v(i,j)]=vr[index_h(i,ord[j])];
                }
                Vr[index_v(iter_max,j)]=0;
            }
            }
        return;
        */
    }
            CHECK(hipMemset(d_G,0,NG*(iter_max+1)*sizeof(float)));
    double relres_tol = relres;
    // return;
    for(outer=0;outer<maxit;outer++){
        initialize(rhs,iter_max+1);initialize(rhsR,iter_max+1);
        //if(outer>=L) kR=0;
        //if(outer>=6) kR=4;
        if(kR>0){
            //GPTLstart ("orth");  
            //
            char dreorth[] = "dreorth";
            orth_mgs(Vr, (iter_max+1), kR, dreorth, P, tid); // reorthogonize Vr -> P_kR 
            orth_vec(P, rR, iter_max+1,  kR, tid); // orthogonize rR to P 
            //GPTLstop ("orth");  
            for(int l=0;l<=iter_max;l++) P[index_v(l,kR)]=rR[l]; //  P_{kR+1} <- [P rR]
            
            swapArrays(d_p, d_G);
            //initialize(h_zp,NG1);
            GPTLstart ("V*P");  
            CHECK(hipMemset(d_p,0,NG*(kR+1)*sizeof(float)));
            CHECK(hipMemcpy(d_P,P,iter_max*(iter_max+1)*sizeof(float),hipMemcpyHostToDevice));
            hipLaunchKernelGGL(gcr_iteration3, numBlocks, threadsPerBlock, 0, 0, iter_max, kR, d_p, d_G, d_P, its, ite, jts, jte, kts, kte, jend);
            GPTLstop ("V*P");  
            GPTLstart ("copy_hp");  
            
            for (is = 0; is < kR; is++) {
                hipblasSdot(h, NG,d_p+kR*NG,1,d_p+is*NG,1, d_xc+is);
            }
            CHECK(hipMemcpy(c1, d_xc, kR*sizeof(float), hipMemcpyDeviceToHost));
            MPI_Allreduce(c1,c2,kR,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
            CHECK(hipMemcpy(d_xc,c2, kR*sizeof(float), hipMemcpyHostToDevice));
            hipLaunchKernelGGL(gcr_iteration3b, numBlocks, threadsPerBlock, 0, 0, iter_max, kR, d_p, d_xc, d_P, its, ite, jts, jte, kts, kte, jend);
            
            hipblasSdot(h, NG,d_p+kR*NG,1,d_p+kR*NG,1, d_xc);
            CHECK(hipMemcpy(c1, d_xc, kR*sizeof(float), hipMemcpyDeviceToHost));
            MPI_Allreduce(c1,c2,1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
            //CHECK(hipMemcpy(d_xc,c2, 1*sizeof(float), hipMemcpyHostToDevice));
            
            hipLaunchKernelGGL(init_p1, num_blocks, block_size, 0, 0, d_p+kR*NG, d_p+kR*NG,(float)1.0/c2[0], its, ite, jts, jte, kts, kte, jend);
            
/*
 #ifdef reorthogize
       
                for (is = 0; is < kR; is++) {
                    tmp_sum=0;
                    para_jend(tmp_sum+=h_p[index4b(i,k,j,kR)] * h_p[index4b(i,k,j,is)];)
                    MPI_Allreduce(&tmp_sum,&beta,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
                    H[index_h(is,kR)] = beta; 
                }
                for (is = 0; is < kR; is++) {
                    para_jend(h_p[index4b(i,k,j,kR)]-=H[index_h(is,kR)]* h_p[index4b(i,k,j,is)];)
                }   

        fprintf(fp,"\nH_kR=\n");
        for(int gi=0;gi<kR;gi++){
            fprintf(fp,"%e ",H[index_h(gi,kR)] );
        }
        fprintf(fp,"\n");
            
                beta= norm2(h_p+kR*NG1, NG1); 
                MPI_Allreduce(&beta,&tmp_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
                tmp_sum=sqrt(tmp_sum);
                para_jend(h_p[index4b(i,k,j,kR)]=h_p[index4b(i,k,j,kR)]/tmp_sum;)
                    
#endif
*/

            //CHECK(hipDeviceSynchronize());
            //CHECK(hipMemcpy(d_G,d_p,NG*(kR+1)*sizeof(float),hipMemcpyDeviceToDevice));
            //hipLaunchKernelGGL(gcr_iteration3b, numBlocks, threadsPerBlock, 0, 0, iter_max, kR, d_G, d_p, d_P, its, ite, jts, jte, kts, kte, jend);
            GPTLstop ("copy_hp");  
            for(is=0;is<=kR;is++){
                for(js=0;js<kR;js++){
                    tmp_sum=0.0;
                    for(int l=0;l<=iter_max;l++){
                    for(int t=0;t<iter_max;t++){
                        tmp_sum+=P[index_v(l,is)]*H0[index_h(l,t)]*P[index_v(t,js)];
                    }}
                    H[index_h(is,js)]=tmp_sum;
                 }
            }
            for(is=0;is<=kR;is++){
                for(js=0;js<kR;js++){
                    H0[index_h(is,js)]=H[index_h(is,js)];
                 }
            }
       
            //GPTLstart ("QR");  
            qr(H, kR+1, kR,  iter_max, Q, R);
            //GPTLstop ("QR");  
            for(is=0;is<=kR;is++){
            for(js=0;js<kR;js++){
                H[index_h(is,js)]=R[index_kr(is,js)];
            }}
                    
            GPTLstart ("V*r0");  
            for(is=0;is<=kR;is++){
                hipblasSdot(h, NG,d_r,1,d_p+is*NG,1, d_xc+is);
            }
            GPTLstop ("V*r0");  
       
            CHECK(hipMemcpy(c2,d_xc, (kR+1)*sizeof(float), hipMemcpyDeviceToHost));
            MPI_Allreduce(c2, rhsR, kR+1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
            for(is=0;is<=kR;is++){
                tmp_sum=0.0;
                for(js=0;js<=kR;js++){
                    tmp_sum+=Q[index_kr(js,is)]*rhsR[js];
                }
                rhs[is]=tmp_sum;
            }
            /*
        fprintf(fp,"\niter=%d,rhs=\n",outer);
        for(int gi=0;gi<=iter_max;gi++){
            fprintf(fp,"%e ",rhs[gi]);
        }
        fprintf(fp,"\n");
        */

        }else{ //rhs(1)=norm(r);V(:,1)=r/rhs(1);
            //beta = norm2(h_r, n);
            hipblasSdot(h, NG,d_r,1,d_r,1, d_xc);
            CHECK(hipMemcpy(c1, d_xc, sizeof(float), hipMemcpyDeviceToHost));
            beta=c1[0];
            #ifdef mpi_reduce
            MPI_Allreduce(&beta,&tmp_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
            beta=sqrt(tmp_sum);
            #else
            beta=sqrt(beta);
            #endif
  
            hipLaunchKernelGGL(init_p1, num_blocks, block_size, 0, 0, d_p, d_r, (float)1.0/beta, its, ite, jts, jte, kts, kte, jend);
            //hipLaunchKernelGGL(init_p2, num_blocks, block_size, 0, 0, d_r0, d_p, its, ite, jts, jte, kts, kte, jend);
            //hipLaunchKernelGGL(init_p1, num_blocks, block_size, 0, 0, d_p, d_r0,(float)1.0/beta, its, ite, jts, jte, kts, kte, jend);
            rhs[0] = beta;rhsR[0] = beta;
        }
    gmres_iter(a0,d_a0, H, H0, d_p, h_zp, d_p1, d_ap, d_ap1, cm_helm, d_cm_helm, rhs, rhsR, js, iter_max, its, ite, jts, jte, kts, kte, jend, tid, kR, Q, ep, flag,  d_r,  buff, d_buff, d_xc, h );
        //    printH0();
        //printH();
        /*
        fprintf(fp,"\niter=%d,rhsR=\n",outer);
        for(int gi=0;gi<=iter_max;gi++){
            fprintf(fp,"%e ",rhsR[gi]);
        }
        fprintf(fp,"\n");
*/
        relres = fabs(rhs[js+1]);/// normb;
        //relres_tol *= relres;

        if (relres  < ep ) {
            flag=0;
        }
        Update(x_c, js, H, iter_max, rhs); //0

        GPTLstart ("update_p");  
        cuinitialize(d_p1,NG1);
        CHECK(hipMemcpy(d_xc,x_c, iter_max*sizeof(float),hipMemcpyHostToDevice));
        hipLaunchKernelGGL(gcr_iteration4b, numBlocks, threadsPerBlock, 0, 0, js, 1, d_xc,  d_p, d_p1 , its, ite, jts, jte, kts, kte, jend);
        GPTLstop ("update_p");  
        #ifdef precondition
        cu_glob_updatehalo_p(h_zp, d_p1, buff, d_buff,its, ite, kts, kte, jts, jte, zero, iter_max,  mytid) ;
        GPTLstart ("precondition");  
        hipLaunchKernelGGL(pacu_svrasr, numBlocks2, threadsPerBlock2, 0, 0, d_p1, 0,d_cm_helm, NX, NY, NZ);
        cu_glob_updatehalo_p(h_zp, d_p1, buff, d_buff,its, ite, kts, kte, jts, jte, zero, iter_max,  mytid) ;
        #endif
        GPTLstop ("precondition");  
        int one=1;
        GPTLstart ("matrixpro");  
#define origin2
#ifdef origin2
//if(kR!=0){
//if(kR==0){
if(outer>=L){ // 5 is best : <5,>=5,>=0,<0
    hipLaunchKernelGGL(cu_matrixpro, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_p1,d_ap, its, ite, jts, jte, kts, kte, jend);
    hipLaunchKernelGGL(gcr_final, numBlocks, threadsPerBlock, 0, 0,  d_x0, d_ap, d_r, d_p1,  its, ite, jts, jte, kts, kte, jend); 
}else{
//#else
    hipLaunchKernelGGL(gcr_final, numBlocks, threadsPerBlock, 0, 0,  d_x0, d_ap, d_r, d_p1, its, ite, jts, jte, kts, kte, jend); 
    //cuinitialize(d_r,NG); //delete
    hipLaunchKernelGGL(cu_matrixpro_x, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_x0,d_r, its, ite, jts, jte, kts, kte, jend, ims, ime, jms, jme);
    hipLaunchKernelGGL(gcr_init, numBlocks, threadsPerBlock, 0, 0, d_f0, d_r, its, ite, jts, jte, kts, kte, jend);
}
#endif
//    hipLaunchKernelGGL(gcr_final1, numBlocks, threadsPerBlock, 0, 0,  d_x0, d_r, d_f0, d_p1, its, ite, jts, jte, kts, kte, jend); 
        GPTLstop ("matrixpro");  
       
        if(kR>0){
            /*
            for(int l=0;l<=iter_max;l++){ 
                tmp_sum=0.0;
                for(i=0;i<iter_max;i++)
                    tmp_sum+=H0[index_h(l,i)]*x_c[i];
                rR[l]=rhsR[l]-tmp_sum;
            }
            fprintf(fp,"\n833rR=\n");
            for(int gi=0;gi<=iter_max;gi++){
                fprintf(fp,"%e ",rR[gi]);
            }*/

 
            harm_Ritz_vec(kR,H0, iter_max, ord, vr,1, tid, rR); 
            /*
        GPTLstart ("dot_p2");  
        CHECK(hipDeviceSynchronize());
        for (is = 0; is <= iter_max; is++) {
            hipblasSdot(h, NG,d_p+is*NG,1, d_r,1, d_xc+is);
        }
            CHECK(hipMemcpy(c1,d_xc, (iter_max+1)*sizeof(float), hipMemcpyDeviceToHost));
        CHECK(hipDeviceSynchronize());
        GPTLstop ("dot_p2");  
            MPI_Allreduce(c1,rR,iter_max+1,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
*/

/*
            fprintf(fp,"\n828rR=\n");
            for(int gi=0;gi<=iter_max;gi++){
                fprintf(fp,"%e ",rR[gi]);
            }
*/  
            for(j=0;j<kR;j++){
                for(i=0;i<iter_max;i++){
                    Vr[index_v(i,j)]=vr[index_h(i,ord[j])];
                }
                Vr[index_v(iter_max,j)]=0;
            }
           /* 
            fprintf(fp,"\nVr=\n");
            for(int gj=0;gj<kR;gj++){
                for(int gi=0;gi<=iter_max;gi++){
                    fprintf(fp,"%3.2E  ",Vr[index_v(gi,gj)]);
                }
                fprintf(fp,"\n");
            }*/

        } 
    hipblasSdot(h, NG,d_r,1,d_r,1, d_xc);
    //hipblasSdot(h, NG,d_ap,1,d_ap,1, d_xc+1);
    //hipblasSdot(h, NG,d_p1,1,d_p1,1, d_xc+1);
    //CHECK(hipMemcpy(h_zp, d_p1, NG1*sizeof(float), hipMemcpyDeviceToHost));
    //CHECK(hipMemcpy(c1, d_xc, sizeof(float), hipMemcpyDeviceToHost));

    hipblasSdot(h, NG,d_p1,1,d_p1,1, d_xc+1);
    CHECK(hipMemcpy(c1, d_xc, 2*sizeof(float), hipMemcpyDeviceToHost));
    //c1[1]=max_norm(h_zp,NG1);
    //beta=c1[0]; //error
    MPI_Allreduce(c1,c2,2,MPI_REAL,MPI_SUM,MPI_COMM_WORLD);
    c2[0]=sqrt(c2[0]);
    c2[1]=sqrt(c2[1]);
    //c2[1]=c2[1];

         //rt_mult_inv(H0, H, Q,  js,iter_max);
        //printQ();
        //com_LOO(Q,Q1,kR,kR);
        //printQ1();
        /*
    CHECK(hipMemcpy(x0, d_x0, NG2*sizeof(double), hipMemcpyDeviceToHost));
    //cuinitialize(d_r,NG); //delete
        matrixpro_x(a0, x0, h_ap,its, ite, jts, jte, kts, kte, jend, ims, ime, jms, jme);
        initialize(h_r,n);
        para_jend(h_r[index3(i,k,j)]  = f0[index3(i,k,j)]- h_ap[index3(i,k,j)];)

    beta = norm2(h_r, n);

    MPI_Allreduce(&beta,&tmp_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    beta=sqrt(tmp_sum);
    */
        //if(tid==12) printf("iter %d,%d RES of gmres is relres=%3.2e,||b-Ax||=%3.2e\n",outer,js,c2[0],beta);
        if(tid==12) printf("iter %d,%d RES of gmres is relres=%3.2e\n",outer,js,c2[0]);
        if(flag==0) break;
        //return;
    }

    CHECK(hipMemcpy(x0, d_x0, NG2*sizeof(double), hipMemcpyDeviceToHost));
    //cuinitialize(d_r,NG); //delete
    /*
    hipLaunchKernelGGL(cu_matrixpro_x, numBlocks_mat, threadsPerBlock_mat, 0, 0, d_a0, d_x0,d_r, its, ite, jts, jte, kts, kte, jend, ims, ime, jms, jme);
    hipLaunchKernelGGL(gcr_init, numBlocks, threadsPerBlock, 0, 0, d_f0, d_r, its, ite, jts, jte, kts, kte, jend);
            hipblasSdot(h, NG,d_r,1,d_r,1, d_xc);
            CHECK(hipMemcpy(c1, d_xc, sizeof(float), hipMemcpyDeviceToHost));
            beta=c1[0];
            #ifdef mpi_reduce
            MPI_Allreduce(&beta,&tmp_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
            beta=sqrt(tmp_sum);
            #else
            beta=sqrt(beta);
            #endif
            */
    //para_jte1(h_zp[index4b(i,k,j,0)]=x0[index_x(i,k,j)];)
    //    int one=1;
//        glob_updatehalo_p_(h_zp, &zero, &one, &its, &ite, &kts, &kte, &jts, &jte);
/*
        matrixpro_x(a0, x0, h_ap,its, ite, jts, jte, kts, kte, jend, ims, ime, jms, jme);
        initialize(h_r,n);
        para_jend(h_r[index3(i,k,j)]  = f0[index3(i,k,j)]- h_ap[index3(i,k,j)];)

    beta = norm2(h_r, n);

    MPI_Allreduce(&beta,&tmp_sum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    beta=sqrt(tmp_sum);
    */

    //CHECK(hipMemcpy(c1, d_xc, sizeof(float), hipMemcpyDeviceToHost));
    //if(tid==5|| tid==12)
    //    printf("iter %d,%d RES of gmres is %e \n",outer,js,beta);

    delete[] ord;
    //delete[] rR;
    //delete[] Vr;
    //delete[] P;
    delete[] G;

    delete[] vr;
    delete[] rhs;
    delete[] rhsR;
    //delete[] H;
    //delete[] H0;
    delete[] x_c;
    //delete[] Q;
    delete[] R;

}

